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Abstract 

When numerically integrating canonical Hamiltonian systems, the long-term 
conservation of some of its invariants, among which the Hamiltonian function 
itself, assumes a central role. The classical approach to this problem has led to 
the definition of symplcctic methods, among which we mention Gauss-Lcgcndrc 
collocation formulae. Indeed, in the continuous setting, energy conservation 
is derived from symplecticity via an infinite number of infinitesimal contact 
transformations. However, this infinite process cannot be directly transferred to 
the discrete setting. By following a different approach, in this paper we describe 
a sequence of methods, sharing the same essential spectrum (and, then, the 
same essential properties), which are energy preserving starting from a certain 
element of the sequence on, i.e., after a finite number of steps. 
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1. Introduction 

In order to make easier the following considerations, it is necessary to stress 
the differences between a continuous problem, for example ([T]), and a discrete 
problem obtained by applying to it a whatever numerical method. The main 
difference between them, very often underrated by many authors, is the lack 
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of continuity in time. The latter does not affect many aspects such as, e.g., 
the study of the qualitative behavior of solutions around asymptotically stable 
critical points (stability analysis). In fact, the respective theories, with minor 
changes, are very similar (see, e.g., [13 )■ As a consequence, many tools, already 
devised in the continuous analysis, can be transferred to the discrete analysis 
almost unchanged. This is the case, for example, of the linearization around 
asymptotically stable critical points, which has been extensively used in the 
numerical analysis of methods for differential problems (linear stability analy- 
sis). There are, however, other aspects for which continuity plays an essential 
role. For example, two solutions of the continuous problem need to stay away 
from each other while, in the discrete case, they may interlace (without having 
common points, of course) . This fact has many mathematical and even physical 



implications (see, e.g., |26[). In this paper we shall deal with another case in 



which continuity plays an essential role. It regards the role of symplecticity 
which is central in discussing energy conservation in continuous Hamiltonian 
problems, while it is less crucial in the energy conservation of discrete problems. 
This depends on the interplay between infinitesimal contact transformations 
and the need of infinite processes (number of iterations) which cannot be opcra- 
tively used in Numerical Analysis. This question has been already discussed in 
previous papers (see, e.g., 0). Here, after a rapid introduction to the subject, 
we shall focus on a particular aspect, although very important, which concerns 
a property of the numerical methods, in the general class of collocation meth- 
ods, which comes out by no more requiring symplecticity but still providing 
conservation of the Hamiltonian functions on a subset of points of the mesh. 
Clearly, this permits to avoid the drift of energy experienced when using many 
numerical methods proposed in the recent literature. 

With this premise, the structure of the paper is the following: in Section[5]we 
recall the basic facts about canonical Hamiltonian problems and the approaches 
used for their numerical solution; one of them, resulting in the recently intro- 
duced class of Hamiltonian Boundary Value Methods (HBVMs), is then sketched 
in Section [3] in Section U we state the main result of this paper, concerning the 
isospectral property of such methods; in Section [5] such property is further gen- 
eralized to study the existing connections between HBVMs and Runge-Kutta 
collocation methods; a few concluding remarks are finally given in Section [6] 



2. Canonical Hamiltonian problems 

Canonical Hamiltonian problems are in the form 

y = JVH(y), y(t ) = y e R 2m , (1) 

where J is a skew-symmetric constant matrix, the Hamiltonian H(y) is assumed 
to be sufficiently differentiable, and the state vector splits into two blocks, y = 
(q T ,p T ) T , q,p € K m where, for mechanical systems, q denote the positions and 
p the (generalized) momenta. Such problems are of great interest in many fields 
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of application, ranging from the macro-scale of celestial mechanics, to the micro- 
scale of molecular dynamics. They have been deeply studied, from the point of 
view of the mathematical analysis, since two centuries. Their numerical solution 
is a more recent field of investigation, where the main difficulty in dealing with 
them numerically stems from the fact that the meaningful isolated critical points 
of such systems are only marginally stable: neighboring solution curves do not 
eventually approach the equilibrium point either in future or in past times. This 
implies that the geometry around them critically depends on perturbations of 
the linear part. Consequently, the use of a linear test equation, which essentially 
captures the geometry of the linear part, whose utility has been enormous in 
settling the dissipative case, cannot be of any utility in the present case. 

It is then natural to look for other properties of Hamiltonian systems that can 
be imposed on the discrete methods in order to make them effective. The first 
property which comes to mind is the symplecticity of the flow tpt '■= Ho ^ y(t) 
associated with ([T]). This property can be described cither in geometric form 
(invariance of areas, volumes, etc.) or in analytical form: 

99tY j(^f£ ]= j_ 



dy ) V d yo 

In one way or the other, it essentially consists in moving infinitesimally on the 
trajectories representing the solutions. Infinitesimally means retaining only the 
linear part of the infinitesimal time displacement 5t. It can be shown that this 
produces new values of the variables q + Sq, p + Sp which leave unchanged the 
value of the Hamiltonian H (q + Sq,p + Sp) = H(q,p) [Infinitesimal Contact 
Transformation (ICT), see [13l . p. 386]). Consequently, since the composition of 
such infinitesimal transformations maintains the invariance, so does an infinite 
number of them. 

It is not surprising that the first numerical attempts to design conservative 
methods have tried to transfer similar arguments to discrete methods, i.e. to 
design symplcctic integrators p3 . [lij (see also the monographs [25 , 21 , 14 1 for 



more details on the subject). A backward error analysis has shown that sym- 
plecticity seems somehow to improve the long-time behavior properties of the 
numerical solutions. Indeed, for a symplectic method of order r, implemented 
with a constant stepsize h, the following estimation reveals how the numerical 
solution y n may depart from the manifold H{y) = -ff(yo) of the phase space, 
which contains the continuous solution itself: 

H(y n )-H(y ) = O(nhe-^h r ), (2) 

where ho > is sufficiently small and h < ho- Relation ([5]) implies that a linear 
drift of the energy with respect to time t = nh may arise. However, due to the 
presence of the exponential, such a drift will not appear as far as nh < : this 
circumstance is often referred to by stating that symplectic methods conserve 
the energy function on exponentially long time intervals (see, for example, [3, 
Theorem 8.1, p. 367]). This is clearly a surrogate of the definition of stability in 
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that the "good behaviour" of the numerical solution is not extended on infinite 
time intervals. 

As matter of fact, since symplecticity requires an infinite sequence of in- 
finitesimal contact transformations it cannot be transferred "sic et simpliciter" 
to the discrete methods, simply because infinite processes are not permitted in 
Numerical Analysis. A more efficient approach would require to design meth- 
ods which avoid the necessity of using infinite processes while preserving the 
constant of motion, i.e., yielding a numerical solution belonging to the manifold 
H{y) = H(yo). In this paper we consider a recently introduced class of methods 
of any high order that provide energy conservation. More specifically, with any 
given order, one can associate an infinite sequences of methods, differing from 
each other for the number of internal mesh points that cover the same time 
window, say [U , U + h] : the more points we include, the better the conservation 
properties of the method. However, we show that, at least for polynomial Hamil- 
tonian functions, such a process of increasing the number of internal points is 
not infinite. In fact we show that there exists a finite value of new added points 
starting from which the method become conservative, whatever is the stcpsizc 
h used. 

The evolution of the approaches to the problem, i.e. to get efficient energy- 
conserving methods, has been slow. As a matter of fact, the first unsuccessful 
attempts to derive energy-preserving Runge-Kutta methods culminated in the 
wrong general feeling that such methods could not even exist [2(| . One of the 
first successful attempts to solve the problem, outside the class of Runge-Kutta 



methods, is represented by discrete gradient methods (see |22| and references 
therein) which are second order accurate. Purely algebraic approaches have 
been also introduced (see, e.g., (HI), without presen ting any energy preserv- 
ing method. A further approach was considered in 23[, where the averaged 
vector field method was proposed and shown to conserve the energy function, 
although it is only second order accurate. As was recently outlined in ap- 
proximating the integral appearing in such method by means of a quadrature 
formula (based upon polynomial interpolation) yields a family of second order 
Runge-Kutta methods. These latter methods represent an instance of energy- 
preserving Runge-Kutta methods f or p olynomial Hamiltonian problems: their 
first appearance may be found in [16j . under the name of s -stage trapezoidal 
methods. Additional examples of fourth and sixth-order conservative Runge- 
Kutta methods (for polynomial Hamiltonians of suitable degree) were presented 
in [13, E^]. AH such energy-conserving methods have been derived by means of 
the new concept of discrete line integral. 

The evolution of this idea eventually led to the definition of Hamiltonian 
Boundary Value Methods (HBVMs) Q, y|, 5|, which is a wide class of methods 
able to preserve, for the discrete solution, polynomial Hamiltonians of arbitrar- 
ily high degree (and then, a practical conservation of any sufficiently differen- 
tiable Hamiltonian). In more details, in Q HBVMs defined at Lobatto nodes 
have been analysed, whereas in Q HBVMs defined at Gauss-Legendre abscis- 
sae have been considered. In the last reference, it has been actually shown 
that both formulae are essentially equivalent to each other, since the order and 



4 



stability properties of the methods turn out to be independent of the abscissae 
distribution, and both methods arc equivalent, when the number of the so called 
silent stages tends to infinity. In this paper this conclusion is further supported, 
since we prove that HBVMs, when cast as Rungc-Kutta methods, arc such that 
the corresponding matrix of the tableau has the nonzero eigenvalues coincident 
with those of the corresponding Gauss-Legendre formula (isospectral property 
of HBVMs). This property will be also used to further analyse the existing 
connections between HBVMs and Runge-Kutta collocation methods. 

3. Hamiltonian Boundary Value Methods 

The arguments in this section are worked out starting from those used in 0, 
13] to introduce and analyse HBVMs. Starting from the canonical Hamiltonian 
problems (TTJ), the key formula which HBVMs rely on, is the line integral and 
the related property of conservative vector fields: 

JJ(yi) - H(yo) = h [ &(t + rhf\7H(<j(t Q + rh))dr, (3) 
Jo 

for any y\ £ K 2m , where a is any smooth function such that 

c(*o)=2/0> a(t + h)=y 1 . (4) 

Here we consider the case where o~(t) is a polynomial of degree s, yielding an 
approximation to the true solution y(t) in the time interval [to, to + h]. The 
numerical approximation for the subsequent time-step, y\, is then defined by 
After introducing a set of s distinct abscissae, 

0< ci,...,c s < 1, (5) 

we set 

Yi = <x(f + Cih), i = l,...,s, (6) 

so that o~(t) may be thought of as an interpolation polynomial, interpolating 
the fundamental stages Yi, i = 1, . . . , s, at the abscissae ([5]). We observe that, 
due to (@|, o-(t) also interpolates the initial condition yo- 

Remark 1. Sometimes, the interpolation at to is explicitly required. In such a 
case, the extra abscissa cq — is formally added to |3p. This is the case, for 
example, of a Lobatto distribution of the abscissae 0/. 

Let us consider the following expansions of a(t) and a(t) for t € [to, to + h]: 

s s PT 

&(t + rh) = ^2 Ti-Pj ( r ). <K*b + rh) =yo + hJ2 7j / Pj 0) dx, (7) 

7=1 J=l J ° 

where {Pj(t)} is a suitable basis of the vector space of polynomials of degree 
at most s — 1 and the (vector) coefficients {7^} arc to be determined. We 
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shall consider an orthonormal polynomial basis on the interval [0, 1] (though, in 
principle, different bases could be also considered @, [H, lj| 17|): 



/ P i {t)P j (t)dt = 6 ij , i,j = l,...,s, (8) 
Jo 

where Sij is the Kronecker symbol, and Pi(t) has degree i — 1. Such a basis can 
be readily obtained as 

Pi(t) = V2i=l Pi-i(t), i = l,...,s, (9) 

with Pi_i{t) the shifted Legendre polynomial, of degree i — 1, on the interval 
[0,1] (see, e.g., [l[). We shall also assume that H(y) is a polynomial, which 
implies that the integrand in is also a polynomial so that the line integral 
can be exactly computed by means of a suitable quadrature formula. It is easy 
to observe that in general, due to the high degree of the integrand function, 
such quadrature formula cannot be solely based upon the available abscissae 
{ci}: one needs to introduce an additional set of abscissae {ci, . . . , c r }, distinct 
from the nodes {c^}, in order to make the quadrature formula exact: 

f &(t Q + Th) T WH(a(t Q + Th))dT = (10) 
Jo 

s r 

A«H*o + Cih) T VH(a(t + ah)) + ^ &<H*o + c*/i) T Vff(a(t + hh)), 

i=l i=X 

where ft, i = 1, . . . , s, and ft, i — 1, . . . , r, denote the weights of the quadrature 
formula defined at the abscissae {a} U {cj}. Then, according to @, we give 
the following definition. 

Definition 1. The method defined by the polynomial o~(t), determined by sub- 
stituting the quantities in ([7J into the right-hand side of (|10j) , and by choosing 
the unknown coefficient {jj} in order that the resulting expression vanishes, is 
called Hamiltonian Boundary Value Method with k steps and degree s, in short 
HBVM(fc,s), where k = s + r. 

According to [lH, the right-hand side of (|10[) is called discrete line integral 
associated with the map defined by the HBVM(fc,s) method, while the vectors 

Yi = a(t + Cih), i = l,...,r, (11) 

are called silent stages: they just serve to increase, as much as one likes, the 
degree of precision of the quadrature formula, but they are not to be regarded 
as unknowns since, from ([7]) and pip , they can be expressed in terms of linear 
combinations of the fundamental stages ((6]). 

Because of the equality ([TO)) , we can apply the procedure described in Def- 
inition [1] directly to the original line integral appearing in the left-hand side. 
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With this premise, by considering the first expansion in ([7}, the conservation 
property reads 



X>J / Pj(T)VH(a(t + rh))dr = 0, 
„•_-, Jo 



3=1 

which, as is easily checked, is satisfied if we impose the following set of orthog- 
onality conditions: 



f Pj(T)JVH(a(t +Th))dT, j = l,..., a. (12) 
Jo 

Then, from the second relation of (|7|) we obtain, by introducing the operator 
L(f; h)a(t + ch) = (13) 

s pC pi 

a(t Q ) + hy2 P,(x)dx P 3 (r)f(<i(t +Th))dT, CG [0,1], 
j =x Jo Jo 

that a is the eigenfunction of L( JV-ff ; h) relative to the eigenvalue A = 1: 

a = L(JVH;h)a. (14) 

According to Q, (|14p is called the Master Functional Equation defining a: it 
characterizes HBVM(fc, s) methods, for all k > s. Indeed, such methods are 
uniquely defined by the polynomial a, of degree s, the number of steps k being 
only required to obtain the exact quadrature formula (1101) . 
To practically compute a, we set (see © and (J7J) 

s 

Yj = a(t + Cjh) = yo + hy^ j ajjjj, i = l,...,s, (15) 

j'=i 

where 

Pj(x)dx, i,j = l,...,s. 



Inserting (jT2"j) into (HH]) yields the final formulae which define the HBVMs class 
based upon the orthonormal basis {Pj}: 



s „i 

Yi = Vo +AV«>j / Pj(T)JVH(a(t + rh)) dr, i = l,...,s. 

3 = 1 J ° 



(16) 



For sake of completeness, we report the nonlinear system associated with 
the HBVM(fc, s) method, in terms of the fundamental stages {Yi} and the silent 
stages {Yi} (see (fTTj) ), by using the notation 

f(y) = JVH(y). (17) 
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It represents the discrete counterpart of ([T5]) , which may be directly retrieved by 
evaluating, for example, the integrals in (|16p by means of the (exact) quadrature 
formula introduced in (|10[) : 

Yi = (18) 

s / s r \ 

= yo + h^an (XlAi , i(ci)/(l r O+X)A^(ai)/(yi) , i = l,.--,«- 

3=1 \i=l i=l / 

From the above discussion it is clear that, in the non-polynomial case, supposing 
to choose the abscissae so that the sums in (|T8"|) converge to an integral as 
r = k — s tends to infinity, the resulting formula is (|16p , which has been named 
oo-HBVM of degree s or HBVM(oo, s) in 0. This implies that HBVMs may be 
as well applied in the non-polynomial case since, in finite precision arithmetic, 
HBVMs are undistinguishablc from their limit formulae (|16p. when a sufficient 
number of silent stages is introduced, so that a practical energy conservation is 
obtained, for k large enough 0, Q, 0, Ell ■ On the other hand, we emphasize 
that, in the non-polynomial case, (|16|) becomes an operative method only after 
that a suitable strategy to approximate the integrals appearing in it is taken into 
account. In the present case, if one discretizes the Master Functional Equation 
(fT3|) - (fT4|) . HBVM(fc,s) are then obtained, essentially by extending the discrete 
problem ([T8"jl also to the silent stages (fTTj) . In more details, by using (fTT)) and 
introducing the following notation: 

{n} = {ci} U {ci}, {u)i} = {ft} U {Pi}, 

yi = o-(t + T i h) 1 fi = f(a(t + nh)), i=l,...,k, 

the discrete problem defining the HBVM(/c, s) method becomes, 

Vi = yo + h>y] / Pj(.x)dx^ u e Pj{Te)fe, i=l,...,k. (19) 

3=1 ^° f=l 

By defining the vectors y = (yf , . . . , y]\) T and e = (1, . . . , 1) T € R fc , and the 
matrices 

n = diag(wi,...,w fc ), l s ,P s £R to (20) 
whose (i, j)th entry are given by 

(D« - / ' Pii?)te, {V s )ij = Pj{n), (21) 



we can cast the set of equations (fT9|) in vector form as 

y = e ® y + h{l a V^n) ® I 2m /(y), 



with an obvious meaning of /(y). Consequently, the method can be regarded 
as a Runge-Kutta method with the following Butcher tableau: 



(22) 



n 














LOl . .. UJk 



In particular, when a Gauss distribution of the abscissae {ti, . . . ,Tk} is con- 
sidered, it can be proved that the resulting HBVM(fc, s) method Q (see also 
Si): 

• has order 2s for all k > s; 

• is symmetric and perfectly A-stable (i.e., its stability region coincides with 
the left- half complex plane, C~ Q); 

• reduces to the Gauss-Legendre method of order 2s, when k = s; 

• exactly preserves polynomial Hamiltonian functions of degree z/, provided 
that 

(23) 



k > 



4. The Isospectral Property 

Wc are now going to prove a further additional result, related to the matrix 
appearing in the Butcher tableau (f2"2"|) . i.e., the matrix 



r,kxk 



k > 



(24) 



whose rank is s. Consequently it has a (k — s)-fold zero eigenvalue. In this 
section, we are going to discuss its essential spectrum, i.e., the location of the 
remaining s nonzero eigenvalues of that matrix. Before that, we state a couple 
of preliminary results: their proofs follow, respectively, from [TBI, Theorem 5.6, 
p. 83] and from the properties of shifted Legendre polynomials (see, e.g., [l[ or 
the Appendix in S). 

Lemma 1. The eigenvalues of the matrix 

( \ -6 

6 o ••• 



with 



V 



'■• -6-1 

6-i o / 



i 



i>i, 



(25) 



(26) 



2V(2j + l)(2j-l) 

coincide with those of the matrix in the Butcher tableau of the Gauss-Legendre 
method of order 2s. 
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Lemma 2. With reference to the matrices in i2U\) - ^2~l]) , one has 



where 





& J 



V 

with the £j defined by i26]) . 

The following result then holds true. 

Theorem 1 (Isospectral Property of HBVMs). For all k > s and for any 

choice of the abscissae {t^} such that the quadrature defined by the weights {wj} 
is exact for polynomials of degree 2s — 1, the nonzero eigenvalues of the matrix 
A in \2J$ coincide with those of matrix i25\) . characterizing the Gauss- Legendre 
method of order 2s. 

Proof For k = s, the abscissae {t{\ have to be the s Gauss-Lcgendre nodes, so 
that HBVM(s, s) reduces to the Gauss Legendre method of order 2s, as already 
outlined at the end of Section [31 When k > s, from the orthonormality of the 
basis, see ([8]), and considering that the quadrature with weights {uJi} is exact 
for polynomials of degree (at least) 2s — 1, one obtains that (see (f20" |) -([2"T T) ') for 
all i = 1 , . . . , s and j = 1 , . . . , s + 1 , 



k „\ 

(VjflV s+1 ) = y2u t Pi(tt)Pj(tt) = / P l (t)P,{t)At 
e=1 Jo 



Sij, 



and, therefore, 



rjnv s+1 = (i s o) 



By taking into account the result of Lemma [21 one then obtains: 

AV S+1 = l s VjCtV s+1 = l s (Is 0) = Vs+iXs (Is 0) = Vs+i ( X s 



-6 
o 



o 



V" 



\ 



Vs+iX s , 



(27) 



J 



with the defined according to ([26]) . Consequently, one obtains that the 
columns of T s +\ constitute a basis of an invariant (right) subspacc of matrix 
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A, so that the eigenvalues of X s are eigenvalues of A. In more detail, the 
eigenvalues of X s are those of X s (see ([25])) and the zero eigenvalue. Then, 
also in this case, the nonzero eigenvalues of A coincide with those of X s , i.e., 
with the eigenvalues of the matrix defining the Gauss-Legcndre method of order 
2s. □ 

It turns out that such methods, in the form here presented (i.e., having 
chosen the polynomial basis (j9])), can be regarded as a generalization of Gauss 
methods, in the sense that, they share the same nonzero spectrum, for all k > 
s. In the limit k — > oo, the same essential spectrum is retained by the limit 
operator (see {331 ). In the case of a polynomial Hamiltonian, such sequence of 
methods starts to be energy-preserving for a finite value of k. Moreover, even 
though for a general Hamiltonian the method becomes energy-preserving in the 
limit k — > oo, nevertheless, when using finite precision arithmetic, the limit 
is practically obtained for a finite value of k, namely as soon as full machine 
precision accuracy is achieved. 

5. HBVMs and Runge-Kutta collocation methods 

By using the previous results and notations, we now further elucidate the 
existing connections between HBVMs and Runge-Kutta collocation methods. 
Our starting point is a generic collocation method with k stages, defined by the 
tableau 



Tl 






A 


Tfc 






UJl ... U! k 



(28) 



where, for i, j = 1, . . . , k: 

A= (aij) = (^j £j(x)dx^j , ujj = J £j(x)dx, 

lj (r) being the jth Lagrange polynomial of degree k — 1 defined on the set of 
abscissae {t,}. Moreover, given a positive integer s < k, and considering the 
matrices defined in (f20|) -(f2~T j) . we consider the matrix 

VsV^Cl G R kxk 

with projects into the s-dimcnsional subspacc spanned by the columns of V s - 
The class of Runge-Kutta methods we are interested in, is then defined by the 
tableau 



Tl 






A = AiVsV^Q) 








UJl UJk 



(29) 
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We note that the Butcher array A has rank which cannot exceed s, because 
it is defined by filtering A by the rank s matrix Vs'PjO,. The following result 
then holds true, which clarifies the existing connections between classical Rungc- 
Kutta collocation methods and HBVMs. 

Theorem 2. Provided that the quadrature formula defined by the weights {oJi} 
is exact for polynomials at least 2s — 1 (i.e., the Runge-Kutta method defined 
by the tableau i29\) satisfies the usual simplifying assumption B(2s)), then the 
tableau V29j) defines a HBVM(k, s) method based at the abscissae {ti}. 

Proof Let us expand the basis {Pi(t), . . . , P s (t)} along the Lagrange basis 
{£j(r)}, j = 1, . . . , k, defined over the nodes Ti, i = 1, . . . , k: 

k 

It follows that, for i = 1, . . . , k and j = 1, . . . , s: 

/ P j (x)dx = y2P j (T r ) i r (x)dx = y2P j (T r )a ir , 

Ja r=l J ° r=l 

that is (sec (EHD (EU and ([2S])). l s = AV S . Consequently, 

AVsTjn = isVjn, 

so that one retrieves the tableau (|2"2")l which defines the method HBVM(/c, s). □ 

The resulting Rungc-Kutta method (|29p is then energy conserving if applied 
to polynomial Hamiltonian systems ([T]), when the degree of H(y) is lower than 
or equal to a quantity, say v, depending on k and s. As an example, when 
a Gaussian distribution of the nodes {r^} is considered, one obtains (|23j) and, 
moreover, HBVM(fc, s) is also related to the Gauss-Legendre method of order 
2k, according to (f2U)) , whose Butcher array coincides with A, with this choice 
of the nodes {n}. 

Remark 2. It seems like the price paid to achieve such conservation property 
consists in the lowering of the order of the new method with respect to the original 
one (|28[) . Actually this is not true, because a fair comparison would be to relate 
method (|22 p - (|29j) to a collocation method constructed on s rather than on k 
stages, since the actual nonlinear system, deriving by the implementation of 
HBVM{k,s), turns out to have dimension s, as it has been shown in Q/. 

Further implications of the isospectral property of HBVMs, among which an 
alternative proof for their order of convergence, may be found in A further 
alternative proof can be found in 0, @| • 
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6. Conclusions 



In this paper, wc have shown that the recently introduced class of energy- 
preserving methods {HBVM(fc, s)}, when recast as Runge-Kutta methods, have 
the matrix of the corresponding Butcher tableau sharing the same nonzero eigen- 
values which, in turn, coincides with those of the matrix of the Butcher tableau 
of the Gauss method of order 2s, for all k > s such that B(2s) holds. 

Moreover, HBVM(fc, s) defined at the Gaussian nodes {ri,...,r/.} on the 
interval [0, 1] are closely related to the Gauss method of order 2k which, although 
symplectic, is not in general energy-preserving. 
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